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Abstract 

The decoupling of a cold relic, during a decaying-particle-dominated cosmological 
evolution is analyzed, the relic density is calculated both numerically and semi-analytically 
and the results are compared with each other. Using plausible values (from the viewpoint 
of supersymmetric models) for the mass and the thermal averaged cross section times 
the velocity of the cold relic, we investigate scenaria of equilibrium or non-equilibrium 
production. In both cases, acceptable results for the dark matter abundance can be 
obtained, by constraining the reheat temperature of the decaying particle, its mass and 
the averaged number of the produced cold relics. The required reheat temperature is, in 
any case, lower than about 20 GeV. 
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1. Introduction 

The enigma of the Cold Dark Matter (CDM) constitution of the universe becomes 
more and more precisely defined, after the recently announced WMAP results [1, 2], which 
determine the CDM abundance, QcBMh 2 , with an unprecedented accuracy: 

ttcvuh 2 = 0.1126l°;S|l (1.1) 

at 95% confidence level. In light of this, the relic density of any CDM candidate \ (i- e > 
which decouples being non relativistic) , Q^h 2 , is to satisfy a very narrow range of values: 

(a) 0.09 < Q^h 2 and (b) n^h 2 < 0.13 (1.2) 

which tightly restricts the parameter space of the theories which support the existence of % 
The most popular of these are the -R-parity conserving supersymmetric (SUSY) theories 
which identify % to the stable lightest SUSY particle (LSP) [3]. According to the standard 
scenario, (i) % decouples from the cosmic fluid during the radiation-dominated (RD) era, (ii) 
being in chemical equilibrium with plasma and (iii) produced through thermal scatterings 
in the plasma. The condition (ii) is satisfied naturally by the lightest neutralino of the 
minimal SUSY standard model (MSSM), which turns out to have the required strength 
of interactions, being weakly interacting. 
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Although quite compelling, this scenario comes across with difficulties, especially when 
it is applied in the context of economical and predictive versions of MSSM. E.g., in most 
of the parameter space of the Constrained MSSM (CMSSM) [4], x turns out to be bino 
and its f^/i 2 exceeds the bound of Eq. (1.2b). Several suppression mechanisms of £ly.h 2 
have been proposed so far: Bino-sleptons [5] and particularly, for large tan/3, bino-stau [6], 
bino-stops [7, 8] or bino-chargino [9] coannihilations and/or A-pole effect [10] can efficiently 
reduce 0,^-h 2 . Also several kinds of non- universality in the Higgs [11] and/or gaugino 
[12, 13] and/or sfermionic sector [14] can help in the same direction, creating additional 
coannihilation effects. As is expected, a more or less tuning of the SUSY parameters 
is needed in these cases, without a simultaneous satisfaction of other phenomenological 
constraints to be always possible (see, e.g. Ref. [2]). On the other hand, Q^h 2 turns out 
to be lower than the bound of Eq. (1.2a) in other models, e.g., in the anomaly mediated 
SUSY breaking model (AMSBM) where \ is mostly wino [15] or in models based on 
SU(5) gaugino non-universality [16], where x can be higgsino. Then we have to invoke 
another CDM candidate [17], in order for the range of Eq. (1.2) to be fulfilled. 

However, this picture can dramatically change, if the standard assumption (i) is lifted. 
Indeed, since there is no direct information for the history of the Cosmos before the epoch 
of nucleosynthesis (i.e., temperatures T > 1 MeV) the decoupling of x can occur not in 
the RD era. As was pointed out a lot years ago [18] and, also, recently [19, 20], the x. 
decoupling can be related to the decay of a massive scalar particle. The modern cosmo- 
particle theories are abundant in such fields, e.g. inflatons [21, 22, 23, 24], dilatons or 
moduli [25, 26], Polonyi field [27], (/-balls [28] (see, also [29]). During their decay, these 
particles perform coherent oscillations "reheating" the universe [30]. This phenomenon is 
not instantaneous [31, 21]. In particular, the maximum temperature during this period 
is much larger than the so-called reheat temperature, which can be better considered as 
the largest temperature of the RD era [19]. Consequently, the "freeze out" of x. could be 
realized before the completion of the reheating. The cosmological evolution during this 
phase is strongly modified as regards the standard one [30], with crucial consequences to 
Q^h 2 calculation [19, 20, 26, 32]. Namely, two types of ^-production emerge, in contrast 
with (ii): The chemically equilibrium (EP) and the non-equilibrium production (non-EP) 
(in both cases, kinetic equilibrium of x's is assumed [19]). In this paper we extend the 
analysis in Ref. [19], lifting also the assumption (iii) of the standard scenario: We include 
the possibility (which, naturally arises even without direct coupling [23, 24]) that the 
decaying particle can decay to x- The problem has already been faced semi-quantitatively 
in Refs [26, 32, 23, 24] and numerically in Ref. [33] for significantly more massive x's. 

Our numerical and semi-analytical analyses are exposed in sees 2 and 3. The obtained 
results are compared with each other in sec. 4. There, we realize, also, a model independent 
application of our findings in the case of SUSY models inspired x masses and cross sections. 
We find that comfortable satisfaction of Eq. (1.2) can be achieved, by constraining the 
reheat temperature to rather low values, the mass of the decaying particle and the averaged 
number of the produced x's, without any tuning of the SUSY parameters. 

Throughout the text and the formulas, brackets are used by applying disjunctive cor- 
respondence and natural units {h = c = ks = 1) are assumed. 
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2. Dynamics of Massive Particle Decay 

We consider a scalar particle 4> with mass m^, which decays with a rate into 
radiation, producing an average number N% of x's with mass m^, rapidly thermalized. 
We, also, let open the possibility (contrary to Ref. [26]) that x's are produced through 
thermal scatterings in the bath. Our theoretical analysis is presented in sec. 2.1 and its 
numerical treatment in sec. 2.2. Useful approximated expressions are derived in sec. 2.3. 

2.1 Relevant Boltzmann Equations 

The energy density of radiation p R and the number densities of 4>, n^, and x, n^, satisfy 
the following Boltzmann equations [21, 26] (we use the shorthand = (m^ — N^m^) jm^): 

+ 3-f/n^ + T^n^ = 0, (2.1) 

p R + 4iTp R - T^A^m^n^ - 2m^{av) (n\ - n^j = 0, (2.2) 

+ 3H njc + {av) (n\ - nf 2 ^j - F^N^ = 0, (2.3) 

where dot stands for derivative with respect to (w.r.t) the cosmic time t, (av) is the 
thermal-averaged cross section of \ particles times velocity and H the Hubble expansion 
parameter which is given by (Mp = 1.22 x 10 19 GeV is the Planck scale): 

H 2 = ^2 (P* + + Px) > with P<t> = ^4> m 4> n <t> and Px = m x n x- ( 2 - 4 ) 
The equilibrium number density of x, ra^ q obeys the Maxwell-Boltzmann statistics: 

nf{x) = j^^mlx^ 2 e' 1/x P 2 (l/x), where x = T/m^, (2.5) 

g = 2 is the number of degrees of freedom of \ an d P n { z ) = 1 + (4n 2 — l)/8z is obtained 
by asymptotically expanding the modified Bessel function of the second kind of order n. 
The temperature, T, and the entropy density, s, can be found using the relations: 

(a) P B =—g P *T 4 and (b) s = — g s * T 3 , (2.6) 

where g p *(T) \g s *(T)] are the energy [entropy] effective number of degrees of freedom at 
temperature T. Their numerical values are evaluated by using the tables included in 
micrOMEGAs [34], originated from the DarkSUSY package [35]. 

Central role to our investigation plays the reheat temperature. Its precise value, T r h, 
can be found, by solving numerically the following [36]: 

P n {T Th ) = P4> (T rh ). (2.7) 

However, following Ref. [19], we prefer to handle reheat temperature as an input parameter, 
Trh, defining it through an analytic formula, which, however approximates fairly (within 
10%) T r h- This can be expressed in terms of T^, using the following [31, 36]: 

l7r 3 g P *{T RH ) Tl u 

~~ 45 Mp' (2 - 8) 

Eq. (2.8) is derived consistently with Eq. (2.16) and differs from the corresponding defini- 
tion in Ref. [19] by a factor 2 (our choice is justified in sec. 4.1). 
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2.2 Numerical Integration 

The numerical integration of Eqs (2.1)— (2.3) is facilitated by absorbing the dilution 
terms. To this end, we convert the time derivatives to derivatives w.r.t the scale factor, R 
[21, 19]. We find it convenient to define the following dimensionless variables: 

U = n<l> R 3 , f R = Pn R 4 and = nfV . (2.9) 

In terms of these variables, Eqs (2.1)-(2.3) become: 



(2.10) 



where 



df<t> _ v f<t> 
dR ~ ^HR' 

§=-<->^ +r ^4 (2 - 12) 

foZ r>-3/2 

H = ] J — ^ I -(A^m (P U + f R /R + mxh) 1/2 . (2.13) 

Since at early times the energy density of the universe is completely dominated by this 
of (p field, m^rij, the system of Eqs (2.10)-(2.12) is solved, imposing the following initial 
conditions: 

f^Ri) = ni Rf and / R (i?j) = = 0, with i? : = m^ 1 . (2.14) 

2.3 Semi- Analytical Approach 

We can obtain a comprehensive and rather accurate approach of the dynamics of the 
decaying-particle-dominated cosmology, following the arguments of Ref. [19]. At the epoch 
before the completion of reheating, T S> TrH) the universe is dominated by the number 
density of <fi which initially has a large value n : . Its evolution is given by the exact solution 
of Eq. (2.1) [30, 31, 36], which at early times, can be written as: 

Inserting this into Eq. (2.4), we obtain: 

H2 = H '{^) ^ With H " = 3^f AwV (2 - 16) 

At early times the last term of the left part side of Eq. (2.11) can be safely ignored and it 
can be easily solved, after substituting Eqs (2.15) and (2.16) in it, with result: 

/r = IfRc RT (R 5/2 - Rl' 2 ) , with / Rc = r A iJ I - 1 m < ^n I . (2.17) 
Combining Eqs (2.17) and (2.6a), we end up with the following relation: 



7T 



'9p* 
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The function T(R/R{) reaches at (see Fig. 1 

— ) =1.48, a maximum value , 



R\ ( 12 \ 1/4 

1.48, a maximum value T max = 0.767 ( — 2 / Rc ) • (2.19) 




For T < T max the second term in the 
second parenthesis in Eq. (2.18) can be elim- 
inated. From the resulting expression, we can 
derive the following useful relations, which "sig- 
nalize" the deviation from the standard RD 
evolution. Namely: 

• Scale factor-temperature relation: 

9% 3 RT^=(^f Rc y 3 Rl . (2.20) 

Recall that in the RD period is the quantity 
gl* 3 (T)RT which remains constant. 

• Expansion rate-temperature relation: 



FIGURE 1: The evolution ofx = T/m^ ver- 
sus R/Ri, derived from the numerical solu- 
tion ofEqs (2.1)-(2.3) (solid lines), from Eq. 
(2.18) (crosses) and from Eq. (2.20) (open 
circles) for the parameters of Fig. 2-(a [b]) 
(normal [light] grey lines, crosses, circles). 



H 



7T 
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9p* Ir c H l T 



(2.21) 



which is obtained by solving Eq. (2.20) w.r.t 
R and inserting it into Eq. (2.16). Note that 
during the RD phase, H is proportional to T 2 
and not to T 4 . 

In Fig. 1, we show x = T/m^ as a function 
of R/Rj for the case of two representative examples which will be analyzed in sec. 4.1. 
The solid lines are obtained by using Eq. (2.6a) and the precise numerical solution of Eqs 
(2.1)-(2.3), whereas crosses [open circles] are from the approximated Eq. (2.18 [2.20]). 
Good agreement is observed for T max > T > Trh in both cases. 



3. Cold Dark Matter Abundance 

The aim of this section is the calculation of Q^h 2 based on the already obtained 
semi-analytical expressions. We assume that %'s are non-relativistic around the 'critical' 
temperatures T* or Tp (see below) and never in chemical equilibrium after the completion 
of reheating (however, see sec. 3.4). The relevant Boltzmann equation is properly re- 
formulated in sec. 3.1. Then two cases are investigated: %'s do (sec. 3.3) or do not reach 
(sec. 3.2) chemical equilibrium with plasma. The condition which discriminates the two 
possibilities is specified in sec. 3.2. 

3.1 Reformulation of the Boltzmann Equation 

Our first step is the re-expression of Eq. (2.3) in terms of the variables Y' cq l = n^ q '/s 8 ^ 3 
in order to absorb the dilution term. Armed with Eqs (2.20) and (2.21), we are able now 



3 Cold Dark Matter Abundance 



6 



to convert the derivatives w.r.t t, to derivatives w.r.t x = T/m^. In fact, supposing that 
g s * (T) and g p * (T) vary very little from their values at Trh and differentiating the constant 
quantity Rs 8 / 9 (see Eq. (2.20) and (2.6b)) w.r.t i, we obtain: 

-8s = 9Hs =► ±- 1 = -8s'/9Hs, (3.1) 

where prime means derivative w.r.t x. Hence, Eq. (2.3) can be rewritten in terms of the 
new variables as: 

Y = -(av) (Y 2 - y cq2 ) s 8 / 3 + T^N^s- ll '\ (3.2) 

Replacing R by x in Eq. (2.15), through Eq. (2.20), we arrive at: 

n^x) = n^x/xj 8 , (3.3) 

with x 1 = Ti/m^ and T\ the temperature corresponding to R\ derived from Eq. (2.20). 
Differentiating Y w.r.t x and replacing x~ l by Eq. (3.1), we obtain: 

8s' 

Y' = Y±- 1 ^Y' = -— mjt Y. (3.4) 

Substituting Eqs (2.6b), (2.21) and (3.2) in Eq. (3.4), this reads: 

Y' = y n (Y 2 -Y^ 2 )x 3 -y N x- 5 , (3.5) 
64 / 2tt^ \ 

where y n (x) = — ( -^rj flf 1 f Rc m\ g a \ (av) (3.6) 

and y N ( x ) = — \JLj H- 1 N^T^ f Rcni x; 8 m^ 12 g ol . (3.7) 
The shorthand g Q \ = g p ^ 2 g% 3 g l J 2 and g Q \ = g p ^ 2 gs} 1 ^ g^J 2 has been used, with [12]: 

9 V 2 = ^=(l + X M- (3.8) 
Y 9p* \ ^9s* J 

3.2 non-Equilibrium Production 

In this case, Y <C Y cq and so, Y 2 — Y cq2 ~ — Y cq2 . Inserting this into Eq. (3.5) and 
integrating from x Y down to x RH = T R n/m^, we arrive at Y(x BH ) = Irh: 

Y RH = Y n (x RU ) + Y N (x BH ), (3.9) 
,2 /0 9\-n/3 



64 a 2 /2vr 2 \~ ' 

with y B (x) = (-) Iff 1 /rc J„(s) m~ 6 (3.10) 

(2/1 / o 2 \ — H/3 

and Yjv(x) = — ( -^_J tff 1 N% J N (x) f Rc m : 12 n t x~ 8 , (3.11) 
where we have defined the quantities (x I can be approximated by infinity, since x l ^> x m ): 

roc roo 

J n (x')= dxx- 10 e- 2 / x P 2 (l/x)g^ 16/3 g ol (av) and J N (x') = dx x~ 5 g ol . (3.12) 

Jx' Jx' 

The maximum x particles production takes place at x* = T^/m^ (=0.212, for constant 
(av)), where the integrand of J n reaches its maximum [19]. Therefore, the condition which 
discriminates the EP from the non-EP of %'s: 

y B (s„) + Y N (x.) - Y cq (x*) | < £ ^° p n - EP (3.13) 



7 



MASSIVE PARTICLE DECAY AND COLD DARK MATTER ABUNDANCE 



3.3 Equilibrium Production 

In this case, we introduce the notion of freeze out temperature, T F = x F m^ [31, 37], 
which assists us to study Eq. (3.5) in the two extreme regimes: 

• At very early times, when x 3> x F , %'s are very close to equilibrium. So, it is more 
convenient to rewrite Eq. (3.5) in terms of the variable A(x) = Y(x) — Y cq (x) as follows: 

A' = -Y cq/ + y n A (A + 2Y cq ) x 3 - y N x~\ (3.14) 

The freeze-out temperature can be defined by 

A(s F ) = 5 F Y cq {x F ) A(x F )(A(x F ) + 2Y e i(x F )) = 5 F (5 F + 2) Y cq2 (x F ), (3.15) 

where 5 F is a constant of order one determined by comparing the exact numerical solution 
of Eq. (3.5) with the approximate under consideration one. Inserting Eqs (3.15) into 
Eq. (3.14) and neglecting A' (since A' <C Y cq/ ), we obtain the following equation, which 
can be solved w.r.t x F iteratively: 

(lny"(x„))' = Vnr S F (S T +2)Y^(x r )x 3 F -y„x; d /Y«'(x r ), with (3.16) 
= ,„,„,(,,) and ( ta y* W )' = ^_»_5£_-^. (3.17) 

• At late times, when x <C x F , Y S> Y cq and so, Y 2 — Y cq2 ~ Y 2 . Inserting this into 
Eq. (3.5), the value of Y at Trh, ^rh = ^(^ RH ) can be found by solving the resulting 
differential equation. We distinguish the cases: 

i. For N% = 0, the integration from x F down to 

^RH 

be made trivially, with 

result: 

Yw=[^r + 4z(¥Pl M|r^m|j F ] , (3.18) 



Yp 15vr V 45 J 
where we have defined the quantities: 



f X F 

J F = dx x 3 g ol (av) and Y F = (5 F + 1) Y cq (x F ), (3.19) 

RH 

/o 2\~8/3 

with Y cq (x) = f^j 5s - 8/3 m- 5 x^ 2 e -V* P 2 (l/ X ). (3.20) 



The choice 5 F = 1.0 T 0.2 provides the best agreement with the precise numerical solution 
of Eq. (3.5), without to cause dramatic instabilities. 

ii. For N% ^ 0, the integration of the resulting equation can be realized numerically. 
However, fixing (av) and g's to their values at x F , an analytic formula can be derived for 
this case, also (note that x RH < x F ): 

Yrh = (-2 + y F tan (tan" 1 2 + y ^ Yp + y p l n \\ (3 21) 

where Yp is given by Eq. (3.20) and the quantity y F = ^—4 — y nF y NF has been defined. 
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3.4 The Current Abundance 

Our final aim is the calculation of the current x relic density, which is based on the 
well known formula [37, 12]: 

= P°x/Pc = m^s Y /p° c , where Yq = n\/s Q , (3.22) 

/3~ = m^n 1 -, is the current % energy density and so [p°] is the entropy [critical energy] 
density today. With a background radiation temperature of T = 2.726 °K, we arrived at 
the final result: 

fl^h 2 = 2.741 x 10 8 Yq m^/GeV. (3.23) 

For the numerical program, lo is determined at a value i?f large enough so as, Yq 
is stabilized to a constant value (with g 7 s fixed to their values at Trh). For the semi- 
analytical calculation, we assume that there is no entropy production for T < Trh- Then, 
it is convenient to single out the cases: 

i. For = 0, there is no x production for x < x RH . Therefore, the present value 
of Y can be estimated reliably at x RH with Yq = Yrh s^^hh) ( m accor d with Ref. [19]). 

ii. For ^ 0, the residual produced (i.e., for x < x RH ) can annihilate 
satisfying an equation analog to Eq. (2.3) but in a RD background, any more, determined 
by Eq. (2.6a) (similarly to the case of Ref. [28]): 

ii x + 3H RD n x + (av) (n\ - nf 2 ^ - T^N^n^ = 0, with H RB = tt^Pr- (3-24) 

The equilibrium distribution n~ q can be safely neglected since, as it turns out, is in any 
case numerically irrelevant. Using the variable Yrd = n x/ s an d the entropy conservation 
law in the RD era, Eq. (3.24) can be rewritten as: 



Y™ = ^- D (^-r^- 2 ). [12] ~f§ £ <3 - 25) 

The post-late times, i.e. for x < x RH , evolution of can be approximated, using the exact 
solution of Eq. (2.1) [30], the entropy conservation law and the time-temperature relation 
[36] in a RD era, as follows (the ratio x BM /x 1 can be safely ignored): 

M x ) = "rh O^rh) 3 gs*{x)/9s*{x BM )e- x ™/ x2 , (3.26) 

where n RH = n^(x RH ) is given by Eq. (2.15). Inserting Eqs (3.26) and (2.6b) into Eq. 
(3.25), this can be cast in the following final form: 

*RD = y w A ~ W e-4 H /*V3 ; ( 3. 27) 
where y BDn (x) = ^J^M P m x gl /2 {av) (3.28) 
and V W n( x ) = \/^(l^) m p n x r n R H \ H 3 ^ 5 gl /2 g^gjJ-(x w ). ( 3 - 29 ) 



Yo can be obtained, by solving numerically Eq. (3.27) from x RH down to 0, with initial 
condition Y R d(x rh ) = Y(x RH ),s 5 / 3 (x RH ), Y(x RH ) being derived from Eq. (3.21) or (3.9). 
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Even in the exceptional case, (see sec. 4.2) where x decouples a little after the completion 
of reheating (i.e., for x < x RH ), Eq. (3.27) can be applied with Yrd(x rh ) = 0, giving 
reliable results. When the first [second] term in the right hand side of Eq. (3.27) evaluated 
at x RH dominates over the second [first] one, an analytical solution of Eq. (3.27) can be 
easily derived, which works frequently well, as we will see in sec. 4.2: 

(a) Yq = jjiD (^rh ) ^) Yq = Yrx>(x ) + ^ RDJV ("^rh ) ^3 30^ 

1 + Urdu ( x rh)^Rd( x rh) x rh ^ 2ex RH J 

where (av) and g's have been fixed at their values at x RH . Practically, since an unavoidable 
discrepancy enters between the input x RH and the output x r h = T^/m^ (i.e., the solution 
of Eq. (2.7)), x RH in Eqs (3.26), (3.27), (3.29) and (3.30) is to be replaced by a matching 
scale x$ = x rh /5rh with <5rh = 0.9 ± 0.3, in order for the solution of Eq. (3.27) to match 
better the solution of the system of Eqs (2.1)-(2.3). Although we did not succeed to achieve 
a general analytical solution of Eq. (3.27), we consider as a significant development the 
derivation of a result for our problem by solving numerically just one equation, instead of 
the whole system above. 

4. Applications 

Our numerical investigation depends on the parameters: 

"if N%, 7r H , m*, (av) . 

Note that the numerical choice of n 1 turns out to be irrelevant for the result of fl^/i 2 . Just 
for definiteness we determine it, through using Eq. (2.16) with H\ = m^, as in the 
simplest model of chaotic inflation [25, 21]. r^, or equivalently Trh, and m<f, can be related 
through a coefficient which includes an effective suppression scale of the interaction of <f>, 
M e g and the coupling of <j) and x, A [26, 32]. However, since we examine the problem from 
cosmological point of view, we prefer to handle and Irh as input parameters without 
any further reference to M e fj and A, which are obviously particle physics model dependent. 
In our scanning we take into account the following bounds: 

(a) 10 3 GeV <m^<8x 10 14 GeV, (b) < 1 and (c) T RH > 0.001 GeV . (4.1) 

The lower bound of Eq. (4.1a) is just conventional, whereas the upper bound is required, 
such as the decay products of the inflaton are thermalized within a Hubble time, through 
2^3 processes [38]. The later is crucial so that Eqs (2.1)— (2.3) are applicable. The bound 
of Eq. (4.1b) comes from the arguments of the appendix of Ref. [26] and this of Eq. (4.1c) 
from the requirement not to spoil the successful predictions of Big Bang nucleosynthesis. 

Model dependent is also the derivation of (av) from and the residual SUSY spec- 
trum. To keep our presentation as general as possible, we decide to treat and {av) 
as unrelated input parameters, choosing plausible (from the viewpoint of SUSY models) 
values for them. Namely, we focus our attention on the ranges: 

(a) 200 GeV <m^< 500 GeV and (b) 10~ 12 GeV~ 2 < (av) < 10" 8 GeV" 2 . (4.2) 



4 Applications 



10 



The lower bound in Eq. (4.2a) arises roughly from the experimental constraints on the 
SUSY spectra (see, e.g., Fig. 23 of Ref. [2]), whereas the upper is imposed in order the 
analyzed range to be possibly detectable in the future experiments (see, e.g. Ref. [39]). On 
the other hand, (av) of the range of Eq. (4.2b) can be naturally produced in the context 
of SUSY models (see, e.g., Fig. 1 of Ref. [20]). The lower value can be derived in the 
case of a bino LSP or sleptonic coannihilations [5] whereas the upper, in the case of a wino 
LSP [28] or coannihilations with squarks or gauginos or of A-pole effect [8]. Note that the 
rc-dependence (which turns out to be important only for EP) in the case of a bino LSP can 
be reliably absorbed, by fixing x in (av) to x F (e.g., if we had posed (av) = 10~ 10 x GeV -2 
for the residual parameters of Figs l-(ai, 82), we would have obtained Q^h 2 ~ 0.57, which 
could, also, be derived by imposing (av) = 10 _10 /21 GeV -2 ). 

In sec. 4.1 we will illustrate the two kinds of x production and in sec. 4.2 we will 
compare the results of our numerical and semi- analytical r^/1 2 calculations. Finally, in 
sec. 4.3 we will present areas compatible with Eq. (1.2). 

4.1 Equilibrium Versus non-Equilibrium Production 

The operation of the two funda- 



Figure 


2-(ai,a 2 ) 


2-(bi,b 2 ) 


Input Parameters 


m m = 10 6 [350 


GeV, (av) = 


10- 10 GeV" 2 


^RH 


5/350 


0.001/350 




1.4 x 10 -7 


1.2 x 10~ 3 


9p* ( X RH ) 


10.88 


3.36 


9s* { X RH ) 


10.85 


3.91 


Output Parameters 


X-Decoupling 


(Y n + Y N )(x*) 


1.4 x 10~ 10 


2.0 x 10~ 18 


Y eq (x*) 


2.8 x 10~ 16 


2.8 x 10~ 16 


Precise Reheating 


^max 


4378.2 


80.4 




4.53/350 


0.0009/350 


P4> (x A ) (GeV 4 ) 


11826.3 


2.3 x 10~ 12 


(R/Rl)rh 


8.06 x 10 14 


1.38 x 10 20 



mental phenomena described in our wo- 
rk (i.e., the x EP [non-EP] and the 
reheating) are instructively displayed 
in Figs 2- (a [b]). Namely, we depict 
in Figs 2-(ai, bi), p l p ] /s (solid [dotted] 
lines) and in Figs 2-(a2,b2), the 4> [ra- 
diation] energy density Pm^a (dashed 
[dot dashed] lines) versus x (we prefer 
P^/s instead of Y,pjjL, since it offers 
a unified description of the evolution 
before and after x RH ). The needed for 
our calculation inputs and some key- 
outputs are listed in Table 1. 

In both cases, we used the same 
771,0 and we fixed and (av) in the 
middle of the ranges of Eq. (4.2), which 
correspond to f^/7 2 = 1-87 with x F = 
22.25/350 for the standard paradigm 
(see sec. 4.3). By adjusting Trh and 
N%, we achieve EP or non-EP, extracting the central value of f^/i 2 in Eq. (1.1). The dis- 
tinguish is realized by explicitly applying the criterion of Eq. (3.13) for = 74.2/350 (see 
Table 1) and is illustrated in Figs 2-(ai,bi). Higher Trh, and consequently (see sec. 4.2), 
lower N% are required for EP. In this case, the solution of Eq. (3.16) is x F = 20.32/350. 
The choice <5rh = 1-15 [0.68] allows us to reach the numerical solution of Eqs (2.1)-(2.3), 
by solving Eq. (3.27) or employing Eq. (3.30b). 

Let us, now, stress on the crucial hierarchies encountered in our examples (see Table 
1) and explain their implications, in conjunction with the analyses of Refs [19, 24, 40]: 



TABLE 1: Input and output parameters for the two 
examples illustrated in Figs 1 and 2-(a, b). 
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FIGURE 2: (ai, bi) The quantity p| /s (solid [dotted] line) and (a 2 , b 2 ) the cj> [radiation] energy 
densities p^y (dashed [dot dashed] line) versus x = T/m^. We take = 10 6 GeV, = 
350 GcV, (av) = 10~ 10 GcV~ 2 and T RH = 5 GeV, N% = 1.4 x 1(T 7 (ai,a 2 ) for EP or T RH = 
0.001 GeV, N% = 1.2 x 10~ 3 (bi, b 2 ) for non-EP. In both cases, we extract tt^h 2 = 0.115. 

• In both cases, T max = x max m^ > T T ^ = x^m^ with x r h being derived by solving 
Eq. (2.7) and corresponding to {R/Rijxh. 3> {R/R\)o (see Fig. 1). Our definition in Eq. 
(2.8) allows us to approach quite successfully x r ^ by x RH (in contrast with Ref. [19]). Our 
T ma x ~ m ]/ 4 Mp 4 T^ , in accord with Ref. [19], turns out to be larger by about one order 
of magnitude than this, estimated in Ref. [24]. 

• In both cases, T max > > Trh/^ f (*)- This means that x 1S first relativistic and 
then becomes non-relativistic, and so, its characterization as cold relic is self-consistent. 

• In Fig 2-(a [b]), we have T max ~ [<]m^. So, our perturbative approach to the <j) 
decay is well defined [24] and possible effects of plasma masses are certainly negligible [40]. 

Let us finally note that hard-hard and hard-soft contributions to VL^h 2 [24] are not 
included in our approach. For the ranges of Eqs (4.1) and (4.2), the hard-hard contributions 
are certainly negligible, whereas the hard-soft are mostly sub-dominant, since 3> m^. 
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4.2 Numerical Versus Semi- Analytical Results 



Fig. 


Ranges of x-Axis Parameters 


^-Pro- 


<5rh 




= 200 GeV 


m ^ = 500 GeV 


duction 






(0.001 - 0.07) GeV 


(0.001 - 0.12) GeV 


non-EP 


0.8 




(0.07 - 2.5) GeV 


(0.12 - 4) GeV 


EP 


1.1 


3-(bi) 


(10 3 - 


10 7 ) GeV 


non-EP 


0.8 


3-(ci) 


lO" 7 


- 10~ 3 


non-EP 


0.8 


3-(a 2 ) 


(0.001 - 0.015) GeV 


(0.0015 - 0.0052) GeV 


EP 


0.75 




(0.4 - 100) GeV 


(1 - 100) GeV 


EP 


0.8 


3-(b 2 ) 


(10 3 - 


10 9 ) GeV 


EP 


0.85 




(10 9 - 10 13 ) GeV 


EP 


1.1 


3-(c 2 ) 


10 -io 


- 10~ 6 


EP 


0.85 




io- 


- 6 -l 


EP 


1.1 



TABLE 2: The type of x-production and the chosen Jrh's for various ranges of the x-axis para- 
meters and 's in Fig. 3. 

The validity of our semi-analytical approach can be tested by comparing its results for 
il^h? with those obtained by the numerical solution of Eqs (2.1)-(2.3). In addition, useful 
conclusions can be inferred for the behavior of ft^h 2 as a function of our free parameters 
and the regions where each %-production mechanism can be activated. 

Our results are presented in Fig. 3. The solid lines are drawn from our numerical code, 
whereas crosses are obtained by employing (with 5-p = 1) Eq. (3.9 [3.21]) for [non] EP and 
solving numerically Eq. (3.27) with a convenient <5rh (comments on the validity of Eq. 
(3.30) are given, too). The type of x-production and the chosen Jrh's for the parameters 
used in Fig. 3 are arranged in Table 2. 

The light [normal] grey lines and crosses correspond to = 200 [500] GeV. In Fig. 
3-(ai,bi,ci [a 2 ,b 2 ,c 2 ]), we fixed (av) = 10~ 12 [10~ 8 ] GeV -2 . We design O^/i 2 versus: 

• T RH , in Fig. 3-(a x [a 2 ]) for iV- = l(r 6 [10~ 3 ] and = 10 6 GeV. Taking into account 
Table 2, also, we deduce that non-EP is accommodated for lower (av) and Trh and higher 
m^, since in these cases, Y n [Y/v] in Eq. (3.10 [3.11]) decreases, facilitating the satisfaction of 
Eq. (3.13) for non-EP. The consideration of N% allows us to produce acceptable results for 
Cl^h 2 even for non-EP, in contrast with Ref. [20], where N% = 0. When ft^h 2 increases with 
Trh, Eq. (3.30b) works well. However, for larger (av) and/or Trh, ^/i 2 decreases, when 
Trh increases (Fig. 3-(a 2 )). None of Eqs (3.30) works in this regime and so, the numerical 
solution of Eq. (3.27) is indispensable. For = 200 [500] GeV, (av) = 10~ 8 GeV" 2 and 
2rh > 8.8 [22.5] GeV, % de-couples for x < x RH (see sec. 3.4). 

• in Fig. 3-(b x [b 2 ]) for N% = 10~ 6 [10~ 3 ] and T RH = 0.05 [5] GeV. In both cases, 
^l^h 2 decreases as increases. Eq. (3.30b [a]) gives reliable results for > 10 9 and 
my. = 500 [200] GeV in Fig. 3-(b 2 ) and for the parameters of Fig. 3-(bi). 

• N%, in Fig. 3-(ci [c 2 ]) for T RH = 0.05 [5] GeV and = 10 6 GeV. We observe 
that £ly.h 2 increases with Eq. (3.30b [a]) gives reliable results for N% < 10~ 6 and 
rriy = 500 [200] GeV in Fig. 3-(c 2 ) and for the parameters of Fig. 3-(ci). 
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FIGURE 3: Q^h 2 versus Trh (ai,a2), (bi,b2) and N% (ci,C2) for fixed (indicated in the 
graphs) N% and m$, and Trh, Trh and correspondingly. We take = 200 [500] GeV 
(light [normal] grey lines and crosses) and (<rv) — 10~ 12 [10~ 8 ] GeV -2 (ai,bi,Ci [a2,b2,C2]). Tie 
solid lines [crosses] are obtained by our numerical code [semi- analytical expressions]. 
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Comparing Fig. 3-(t>i [b2]) and Fig. 3-(ci [02]), it can be induced that Q^h 2 remains 
invariant for constant N^m^ 1 . This can be understood by the observation that y N in Eq. 
(3.7) is proportional to this quantity. The adjustment of <5rh turns out to be crucial, in 
order for the results from the semi- analytical treatment to approach the numerical ones. 
However, this kind of uncertainties is unavoidable to such approximations (see Ref. [19]). 
Note that the accuracy of our approximations increases as N^m^ 1 is reduced. 

4.3 Allowed Regions 

Requiring Q,^h 2 to be confined in the cosmologically allowed range of Eq. (1.2), one 
can restrict the free parameters. The data is derived exclusively by the numerical program. 
Our results are presented in Fig. 4. The light [normal] grey regions are constructed for 
= 200 [500] GeV. In Fig. 4-(ai,bi,ci [a 2 ,b 2 ,c 2 ]), we fixed (av) = 10~ 12 [10~ 8 ] GeV~ 2 . 
We display the allowed regions on the: 

• Ny. — Trh plane, in Fig. 4-(ai,a2) for = 10 6 GeV. Since Q^h 2 increases with 
Trh for low (av) 's, the increase of Trh requires diminution of N% in Fig. 4-(ai). On 
the contrary, in Fig. 4-(a2), when decrease of Q^h 2 is achieved with increase of Trh for 
Trh ^ 2 GeV, augmentation of N% is needed for increasing Trh- 

• Ny. — plane, in Fig. 4-(bi [b2]) for Trh = 0.05 [5] GeV. Since in both cases, VL^h 2 
in [de] -creases as [m^] increases, increase of entails increase of m^. 

• T RH - plane, in Fig. 4-(ci [c 2 ]) for N% = 10~ 6 [10~ 3 ]. Since for T RH < 2 GeV, 
^l^h 2 increases with Trh (see Fig. 3-(ai, 32)) m<i> is to increase with Trh, as in Fig. 4-(ci). 
Similar region exists, also, for the parameters used in Fig. 4- (02). However, for the sake 
of illustration, in Fig. 4-(c2) we concentrated on the regions with Trh ^ 2 GeV. In these, 
increase of m<^ is dictated with decrease of Trh> since the latter causes increase of fi^/i 2 , 
as shown in Fig. 3-(a2). The upper [lower] bounds of the allowed areas are derived from 
the saturation of Eq. (1.2a [b]) (inversely to all other cases). 

Note, finally, that in the standard paradigm, VL^h 2 is (m^, N%, TRH)-independent and 
turns out to be almost 143 - 146 [0.022 - 0.023] for (av) = W' 12 [10^ 8 ] GeV~ 2 and m LS p 
varying within the range of Eq. (4.2a). This is easily extracted by solving Eq. (2.3) with 
N% = and only a RD background in Eq. (2.4) with g p * fixed at x F m^ with 17 < x^ 1 < 26. 
We checked that this result agrees with this obtained by solving Eq. (76) of Ref. [12]. 

5. Conclusions-Open Issues 

We considered a deviation from the standard cosmological situation according to which 
the CDM candidate, \ decouples from the plasma (i) during the RD era (i.e. after reheat- 
ing) (ii) being in equilibrium (iii) produced through thermal scatterings. On the contrary, 
we assumed that x decoupling occurs (i') during a decaying-massive-particle, (ft, domi- 
nated era (and mainly before reheating) (ii') being or not in chemical equilibrium with the 
thermal bath (iii') produced by thermal scatterings and directly from the <j> decay. 

We solved the problem (i) numerically, integrating the relevant system of the differen- 
tial equations (ii) semi-analytically, producing approximate relations for the cosmological 
evolution before reheating and solving the properly re-formulated Boltzmann equations. 
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FIGURE 4: Regions allowed byEq. (1.2) on the Trh-A^ plane (ai, a 2 ) for = 10 6 GeV, m^-N^ 
plane (bi [b 2 ]) forT^ = 0.05 [5] GeV and m -T RH plane (ci [c 2 ]) for = 1Q- 6 ^. Wetakem^ = 
200 [500] GeV (light grey [normal grey] areas) and (crv) = 10~ 12 [10~ 8 ] GeV" 2 (ai,bi,Ci [a 2 ,b 2 ,c 2 ]). 
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The second way facilitates the understanding of the problem and gives, in most cases, 
sufficiently accurate results, provided a suitable <5rh is chosen. The variation of S^/i 2 w.r.t 
our free parameters (m^, N%, Trh) was investigated and regions consistent with the present 
CDM bounds are constructed, using m^'s and {av)'s commonly allowed in SUSY models. 

These scenaria obviously let intact the SUSY parameter space but require rather low 
2rh- This can be accommodated in AMSBM [26], in models with intermediate scale 
unification [32] and within the context of (/-balls decay [28]. Also, low Trh naturally arises in 
models of thermal inflation [41] which, as a bonus, may overcome the problem of unwanted 
relics (e.g., gravitino, moduli). Finally, restrictions on Trh arising from baryogenesis and 
neutrino cosmology have been, also, studied in Refs [42, 43]. 

Our formalism can be easily extended to include coannihilations and pole effects. 
Therefore, it can become applicable for the calculation of Q^h 2 in the context of specific 
SUSY models. Also, these scenaria can assist us to the reduction of h? in cases, where 
it turns out to be even more enhanced than in the standard scenario, as in the presence 
of Quintessence [44]. Similar analysis of the ^-decoupling during the extra dimensional 
cosmological evolution [45] may be also, possible. 
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